%%%% equilibrium outcomes with kink solution

function [K_1_grid, tau_k_2_grid, L_1_grid, L_2_grid, q_1_grid, B_0_grid, tau_l_2_grid] = eqm_kink(B_e_0, par, grid_number, PSI_step, PSI_start)

K_1_grid     = zeros(grid_number,1); 
L_1_grid     = zeros(grid_number,1); 
L_2_grid     = zeros(grid_number,1); 
tau_k_2_grid   = zeros(grid_number,1);
tau_l_2_grid   = zeros(grid_number,1);
q_1_grid     = zeros(grid_number,1);
B_0_grid     = zeros(grid_number,1);

K_star = B_e_0 /(1 - par.phi);

for jj = 1:1:grid_number
    PSI = PSI_step * (jj - 1) + PSI_start;
    
    K_1 = K_star;

    L_1 = (1 + PSI) / (1 + (1 + par.NU) * PSI) * par.A / par.MU;
    
    L_2 = (par.A * (1 - par.ALPHA) / par.MU * (1 + PSI) / (1 + (1 + par.NU) * PSI))^(1 / (par.ALPHA + par.NU))...
        * K_1^(par.ALPHA / (par.ALPHA + par.NU));
    
    q_1 = 1;
    r_2 = par.ALPHA * par.A * (K_1 / L_2)^(par.ALPHA - 1) + 1 - par.DELTA;
    tau_k_2 = 1 - q_1 / par.BETA / r_2;
    w = par.A * (1 - par.ALPHA) * (K_1 / L_2)^par.ALPHA;
    tau_l_2 = 1 - L_2^par.NU * par.MU / w;
    
    C_1     = par.A * L_1 - K_1; % goods mkt clearing in period 1
    C_2     = par.A * K_1^par.ALPHA * L_2^(1 - par.ALPHA);
    B_0     = (C_1 - par.MU * L_1^par.NU * L_1) + par.BETA * (C_2 - par.MU * L_2^par.NU * L_2) - (1 / par.phi - 1) * max(0, K_1 - B_e_0 / (1 - par.phi));


    K_1_grid(jj) = K_1;
    L_1_grid(jj) = L_1;
    L_2_grid(jj) = L_2;
    tau_k_2_grid(jj) = tau_k_2;
    tau_l_2_grid(jj) = tau_l_2;
    q_1_grid(jj) = q_1;
    B_0_grid(jj) = B_0;
end